   function [V_w,A_polw,W_polw] = vworkz_tsp(A,W,last,a,y,lastfuncW,solver_param,f_N,mu,rho,sigma_v,phi,gamma,taxfz,mortp_80,z)
   
    sigma  = solver_param(1);
    beta   = solver_param(2);
    aT     = solver_param(3);
    aS     = solver_param(4);
    r      = solver_param(5);
    V_term = solver_param(6);
    scaler = solver_param(7);
    mortp  = mortp_80(aT-aS-a);
    itau   = taxfz{z,1}; 
    ptau   = taxfz{z,2};
    tau_k  = 0 ;
	
	%Inputs:
	%A = state space for assets age 50
	%W = state space for DC balance age 50
	%f = state space for labor disutility
	%last = 1 if computing last period functions
	%ptau = payroll tax rate 
    %itau = income tax rate
    %tau_k = tax on capital     
	%r = interest rate
	%aT = Terminal age
	%aS = start age
	%a = current age
	%y = income path
	%lastfuncW/R = previously estimated value and policy functions
    %i_z = index of freeze age 
	%V_term = terminal value (=0)
    %scaler = consumption utility scaler
	%f_N = dimension of f state space for rouwenhorst approximation
	%phi = labor disutlity slope on age term
	%gamma = labor disutility constant term
	%rho = labor disutility persistence
	%sigma_v = labor disutility idiosyncratic variance
	
	%Outputs 
	%V_w = value function when working 
	%A_polw = policy function for saving non-DC assets
	%W_polw = policy function for saving DC assets
   
   
	%Parameters for disutility grid (Rouwhenhorst approximation)
	[f,P] = rouwenhorst(f_N,mu,rho,sigma_v);

	% Pre-create state space grid and empty value and policy functions 
    W_0 = W.*(1+r)^(aT-aS-a-1); % DC balance grid current period
    W_1 = W.*(1+r)^(aT-aS-a);   % DC balance grid next period
    A_0 = A.*(1+r)^(aT-aS-a-1); % Non-DC grid current period
    A_1 = A.*(1+r)^(aT-aS-a);   % Non-DC grid next period
	[A_nd,W_nd] = ndgrid(A_1,W_1);  
	V_w = zeros(length(A_1),length(W_1),length(f));
	A_polw = zeros(size(V_w));
	W_polw = zeros(size(V_w));
	
	
	%Solve for policy functions and value function
	for i = 1:length(A_0) %loop over current A space 
	    for j = 1:length(W_0) %loop over current W space
            for k = 1:length(f) %loop over disutility of labor space 
            
                %Set income,retirement annuity, and continuation values 
                if last == 0
                    yt = y(aT-aS-a);
                    itau_l = itau((aT-aS-a-1)*(length(W_0)^2)+1:(aT-aS-a)*length(W_0)^2);
                    itau_l = reshape(itau_l,[length(W_0),length(W)]);
                    ptau_l = ptau((aT-aS-a-1)*(length(W_0)^2)+1:(aT-aS-a)*length(W_0)^2);
                    ptau_l = reshape(ptau_l,[length(W_0),length(W_0)]);  
                    %V_prev = V from prior iteration/future age (i.e. age+1)
                    V_prev = lastfuncW;	
                    %Compute continuation value expectation term 
                    EV = zeros(size(A_nd));
                    for l = 1:length(f)
                        temp = V_prev(:,:,l).*P(k,l);
                        EV = EV+temp;
                    end
                    V_cont1 = EV;
                elseif last == 1
                    yt = y(end);
                    itau_l = itau((aT-aS-1)*(length(W_0)^2)+1:(aT-aS)*length(W_0)^2);
                    itau_l = reshape(itau_l,[length(W_0),length(W_0)]);
                    ptau_l = ptau((aT-aS-1)*(length(W_0)^2)+1:(aT-aS)*length(W_0)^2);
                    ptau_l = reshape(ptau_l,[length(W_0),length(W_0)]);   
                    V_cont1 = V_term;
                end

                %Determine consumption (there are three states for DC
                %matching)            
                %first deduct taxes
                c_hold1 = zeros(size(W_nd));
                c_hold2 = zeros(size(W_nd));
                c_hold3 = zeros(size(W_nd));
                inc1 = yt.*(1-((W_nd./(1+r)-W_0(j))./(2.*yt)-(.01/2))); %income less DC contrib type 1 
                inc2 = yt.*(1-((W_nd./(1+r)-W_0(j))./(1.5.*yt)-(.01/1.5))); %income less DC contrib type 2 
                inc3 = yt.*(1-((W_nd./(1+r)-W_0(j))./(yt)-.05)); %income less DC contrib type 3             
                itax = (1-itau_l(j,:)'); %income tax rate
                ptax = yt.*ptau_l(j,:)'; %payroll tax liability
                for l = 1:length(W)
                    c_hold1(:,l)=inc1(:,l).*itax(l)-ptax(l);
                    c_hold2(:,l)=inc2(:,l).*itax(l)-ptax(l);
                    c_hold3(:,l)=inc3(:,l).*itax(l)-ptax(l);
                end
                %Add in non-pension saving
                cons1 = c_hold1+A_0(i)-A_nd./(1+r*(1-tau_k));
                cons2 = c_hold2+A_0(i)-A_nd./(1+r*(1-tau_k));
                cons3 = c_hold3+A_0(i)-A_nd./(1+r*(1-tau_k));
                T=(W_nd./(1+r)-W_0(j))./yt;
                T1= (T<.07);
                T2= (T<.10 & T>=.07);
                cons = cons1.*T1 + cons2.*(T2) + cons3.*(1-T1-T2); 
                cons(cons<0)=0; %zero out negative consumption levels

                %Restrict W accumulation to 54500 (2010 tax rules)
                T=(W_nd./(1+r)-W_0(j)<=54500); 
                cons = cons.*T; 

                %Value function
                d = (phi*(aT-a)+gamma+f(k));
                    v_w_temp = scaler.*(1/(1-(1/sigma))).*(cons).^(1-(1/sigma))-d...
                      +(1-mortp).*beta.*V_cont1 + 500;
                %Disqualify drawing down W when working
                v_w_temp(:,1:j-1)=-Inf;
                %Eliminate neg consumption levels from maxima
                T = (cons>0);
                v_w_temp = T.*v_w_temp + (1-T).*-1e50;

                %Determine maxima 
                %max over W dimension
                [dimwmax,ind1]=max(v_w_temp,[],2);
                %max over A dimension
                [vmax,ind2]=max(dimwmax);
                %index of W policy (pick ind2nd row from ind1st column)
                Wind = ind1(ind2);       
                %index of A policy
                Aind = ind2;

                %assign value and policy functions
                V_w(i,j,k)=vmax; 
                A_polw(i,j,k)=A_1(Aind);
                W_polw(i,j,k)=W_1(Wind);
            end
	    end
	end
	
   
   end
